library(sf)
library(tidyverse)

#インタラクティブな地図作成
library(mapview)

目的

ある地物（地上にある物）や地点周辺の統計地図の作成，および統計量（合計値，平均値など）の計算．ここでは，多摩ニュータウン（ポイントデータ）周辺の人口（メッシュデータ）を可視化．また多摩ニュータウンから500mの円を描き，その円と重なるメッシュデータの人口数を計算．

知りたい統計量の例

多摩ニュータウンと多摩市周辺の人口分布

ダウンロード

行政区域

  • 国土数値情報ダウンロードサービス　>　行政区域（ポリゴン）　>　東京都　>　世界測地系　令和4年
    • ダウンロードしたZIPファイルは同じフォルダ（ディレクトリ）に納め，ZIPファイルを選択 >　右クリック　>　すべて展開　>　展開先の選択とファイルの選択　>　展開を実行する．展開後は最初にダウンロードしたZIPファイルを削除．
    • シェープファイルの読込（read_sf()）：行政区域をTokyo_mapと名付ける．
      • 多摩市周辺の行政区域の抜き出し（filter()）：Tama_shi_mapと名付ける．
#東京の行政区域（令和4年）
Tokyo_map<-
  read_sf("N03-20220101_13_GML/N03-22_13_220101.shp")

#多摩市周辺（多摩，町田，日野，稲城，八王子）
Tokyo_map %>%
  filter(N03_007=="13224" | N03_007=="13209" |
           N03_007=="13212" | N03_007=="13225" |
           N03_007=="13201") ->
  Tama_shi_map

#多摩市の行政区域の可視化
ggplot()+
  geom_sf(data=Tama_shi_map)

ニュータウン

  • 国土数値情報ダウンロードサービス　>　ニュータウン（ポイント）　>　東京都　>　 世界測地系　平成25年
    • シェープファイルの読込（read_sf()）：行政区域をTokyo_newtownと名付ける．
    • missing crsのエラーが出るため，crs（座標参照系）を設定．
      • WGS84：1984年に定められた世界測地系．
    • 多摩ニュータウンを抜き出し，Tama_newtownと名付ける．
#東京ニュータウン
Tokyo_newtown<-
  read_sf("P26-13_13/P26-13_13/P26-13_13.shp", crs="WGS84")

#多摩ニュータウンの選抜
Tokyo_newtown %>%
  filter(P26_005=="多摩ニュータウン") ->
  Tama_newtown

統計データ：500mメッシュ

市区町村別メッシュ・コード一覧（残念ながら3次メッシュの情報だけの模様）から多摩市のコードが5339からはじまることがわかる．

人口

  • 統計地理情報システム　>　統計データダウンロード　>　国勢調査　>　2020年　>　4次メッシュ（500mメッシュ）　>　人口及び世帯　>　都道府県で絞込みはコチラ　>　13 東京都　>　M5339 1　
  • read.table()を用いてデータを読込．Tokyo_popと名付ける．
#人口
Tokyo_pop<-read.table("tblT001101H5339/tblT001101H5339.txt",
                      header=TRUE, sep=",")

#便宜上1行目削除（データ名が2行にわたるため）
Tokyo_pop %>% 
  slice(-1) ->
  Tokyo_pop

#メッシュデータと結合するため，KEY_CODEのデータ型変換
Tokyo_pop %>%
  mutate(KEY_CODE=as.character(KEY_CODE)) ->
  Tokyo_pop

境界

メッシュデータの人口総数を地図上に可視するには，シェープファイルと結合する必要がある．そこで，メッシュデータの境界データをダウンロードし，結合する．

  • 統計地理情報システム　>　境界データダウンロード　>　4次メッシュ（500mメッシュ）　>　世界測地系緯度経度・Shapefile　>　都道府県で絞込みはコチラ　>　13 東京都　>　M5339

    • read.sf()を用いてデータを読込．meshと名付ける
    • meshとTokyo_popをKEY_CODEで結合．
mesh<-
  read_sf("HDDSWH5339/MESH05339.shp")

#データ結合
Tokyo_mesh<-
  left_join(mesh, Tokyo_pop, by=c("KEY_CODE"))

#人口総数（T001101001）を数値（integer：整数）に変更
Tokyo_mesh %>%
  mutate(T001101001=as.integer(T001101001)) ->
  Tokyo_mesh

地図上に可視化

Tokyo_meshは東京都全体の地図．ここでは，多摩市周辺だけを示したい．そこで，多摩市周辺の行政区域（Tama_shi_map）とTokyo_meshの重複する部分を抽出する．

  • st_intersection()を用いて，2つのシェープファイル（Tokyo_meshとTama_shi_map）の重複部分を抽出．
    • そのためには以下の準備で示されている指示が必要．2
#準備
proj<-
  "+proj=longlat +datum=WGS84"

Tokyo_mesh_proj<-
  Tokyo_mesh %>%
  st_transform(proj) 

Tama_shi_proj<-
  Tama_shi_map %>% 
  st_transform(proj) 

st_agr(Tokyo_mesh_proj)="constant"
st_agr(Tama_shi_proj)="constant"

#多摩周辺の地図に人口が加わったメッシュデータの完成
Tama_mesh<-
  st_intersection(x=Tokyo_mesh_proj, y=Tama_shi_proj)

多摩ニュータウンと多摩周辺の人口分布の可視化

ggplot()

  • NA（欠損値）は情報（ここでは人口総数）が含まれないことを意味する．
#完成図
ggplot()+
  geom_sf(data=Tama_mesh, aes(fill=T001101001))+
  scale_fill_viridis_c(option="G", direction=-1)+
  geom_sf(data=Tama_newtown, fill="yellow", color="yellow")+
  labs(fill="人口総数",
       caption="出典：国土交通省国土数値情報")+
  ggtitle("多摩NT位置と多摩市周辺人口（2020年）")+
  theme_bw()

mapview()

mapview()を用いる利点：マウスカーソルを自ら動かすことで知りたいメッシュ内の人口総数を調べられるようになる．3

  • layer.nameを用いて凡例を変更．
#viridisカラーの使用（好み）
library(viridis)

#色の指示（好み）
my_color<- 
  viridis(7, option="mako", direction=-1)

#mapview表記
#col.regions=my_colorは好み（削除しても問題ない）
mapview(Tama_mesh, zcol="T001101001", 
        alpha.regions=0.6, col.regions=my_color,
        layer.name="人口総数")+
  mapview(Tama_newtown, cex=2.5,
          color="yellow", col.regions="yellow")

周辺500mの可視化

#単位を付与（ここでは周辺500m）
library(units)

Tama_newtown %>%
  st_buffer(geometry, dist=set_units(500, m)) ->
  buffer

#可視化1（500mの円）
ggplot()+
  geom_sf(data=Tama_mesh, aes(fill=T001101001))+
  scale_fill_viridis_c(option="G", direction=-1)+
  geom_sf(data=buffer, color="red", fill="NA")+
  geom_sf(data=Tama_newtown, fill="yellow", color="yellow")+
  labs(fill="人口総数",
       caption="出典：国土交通省国土数値情報")+
  ggtitle("多摩ニュータウン周辺500mの範囲")+
  theme_bw()

人口総数（の合計値）の計算

Buffer_Pop<-
  lengths(st_intersects(Tama_mesh, buffer))>0 
Tama_mesh %>%
  filter(., Buffer_Pop==TRUE) ->
  Tama_mesh2
#可視化2（円を重複するメッシュ（赤）に変更）
ggplot()+
  geom_sf(data=Tama_mesh, aes(fill=T001101001))+
  scale_fill_viridis_c(option="G", direction=-1)+
  geom_sf(data=Tama_mesh2, color="red", fill="NA")+
  geom_sf(data=Tama_newtown, fill="yellow", color="yellow")+
  labs(fill="人口総数",
       caption="出典：国土交通省国土数値情報")+
  ggtitle("多摩ニュータウン周辺メッシュ")+
  theme_bw()

周辺500mのメッシュに含まれる人口総数の計算

#重複するメッシュ（KEY_CODE）を取り除く
#（.keep_all=TRUE）その他の変数残す
Tama_mesh2 %>% 
  distinct(KEY_CODE,.keep_all=TRUE)->
  Tama_mesh3

#欠損値無視オプション（na.rm）
Tama_mesh3 %>%
  summarize(sum_pop=sum(T001101001, na.rm=TRUE))
sum_pop geometry
147078 MULTIPOLYGON (((139.375 35…. 
#上記より下記のよう示す方がよいかも
sum(Tama_mesh3$T001101001, na.rm=TRUE)
## [1] 147078

mapview()で確認

mapview()を用いて，周辺500mと重なるメッシュのみを可視化．

#mapview
mapview(buffer, color="red",
        alpha.regions=1, lwd=1.2,
        legend=FALSE, homebutton=FALSE)+
  mapview(Tama_mesh2, zcol="T001101001", 
          alpha.regions=0.6, col.regions=my_color,
          layer.name="人口総数", homebutton=FALSE)+
  mapview(Tama_newtown, cex=2.5,
          color="yellow", col.regions="yellow",
          legend=FALSE, homebutton=FALSE)

謝辞

ポリゴンと重なるメッシュの検索，抽出方法について山本雅資氏にRコードを教示いただいた．

Rによる地理空間データの可視化

  1. 統計地理情報システムの国勢調査からは1995年から2020年（5年ごと）までの人口総数（男女），世帯総数などがわかる（アクセス日：2022-10-29）．2015年からは年齢別人口，外国人人口など詳しく記録されている．↩︎

  2. 準備に示されているコードはコミュニティバス路線上の人口密度の可視化に説明あり．↩︎

  3. mapview()を用いた地理空間データの可視化についてはインタラクティブな地図によるデータの可視化参照．↩︎